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Abstract 

We propose an alternative interpretation of Markovian transport models based 
on the well-mixed condition, in terms of the properties of a random velocity field 
with second order structure functions scaling linearly in the space time increments. 
This interpretation allows direct association of the drift and noise terms entering 
the model, with the geometry of the turbulent fluctuations. In particular, the well 
known non-uniqueness problem in the well-mixed approach is solved in terms of 
the antisymmetric part of the velocity correlations; its relation with the presence of 
non-zero mean helicity and other geometrical properties of the flow is elucidated. 
The well-mixed condition appears to be a special case of the relation between con- 
ditional velocity increments of the random field and the one-point Eulerian velocity 
distribution, allowing generalization of the approach to the transport of non-tracer 
quantities. Application to solid particle transport leads to a model satisfying, in the 
homogeneous isotropic turbulence case, all the conditions on the behaviour of the 
correlation times for the fluid velocity sampled by the particles. In particular, cor- 
relation times in the gravity and in the inertia dominated case, respectively, longer 
and shorter than in the passive tracer case; in the gravity dominated case, correla- 
tion times longer for velocity components along gravity, than for the perpendicular 
ones. The model produces, in channel flow geometry, particle deposition rates in 
agreement with experiments. 
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I. Introduction 



In Lagrangian models, the concentration of a quantity transported by a turbulent 
field is reconstructed from the trajectories of the individual particles advected by the 
flow m 12 El El- Since each trajectory is associated with an independent realization 
of the turbulent flow, what is obtained is actually the mean concentration profile, 
given a certain distribution of sinks and sources for the transported quantity. 

One keeps into account the fact that the turbulent length and time scales are 
non-zero, by assuming that the particle velocity obeys a Langevin equation; this 
results in the system of equations, in terms of the particle Lagrangian velocity v 
and coordinate x: 

dx = vdt 

(iv = a(v, X, t)dt + dw (1) 
{dw'dw^) = S'^(v,x,t)dt 

A strong physical motivation for a Lagrangian model in the form of Eqn. (^J is the 
linear scaling of the Lagrangian velocity structure functions for inertial range time 
separations [5J. In the case of passive tracers, we have 

{[v\t) - v\0)][v\t) - v^{0)]) ~ S'^Coet (2) 

with e the mean viscous dissipation and Cq a universal constant; Eqn. will then 
result from assuming a Markovian behaviour for v, i.e. that a* and B^^ depend solely 
on the current values of v and x and not on their previous history. 

The well-mixed condition, introduced in lead to a great advance in La- 
grangian models, providing a simple technique for expressing the drift coefficient 
a in Eqn. in terms of observed properties of the flow. 

In the current approach, it is assumed that the noise term dw in Eqn. 
accurately represents, in high Reynolds number turbulent regimes, the inertial range 
scaling of the Lagrangian velocity increments: 

F^(v,x,t) = 5^^Coe(x,t) (3) 

Then, at least when the flow is incompressible, the well-mixed criterion allows us to 
determine the drift coefficient a in terms of e and the Eulerian Probability Density 
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Function (PDF) for the fluid velocity u(x, t): ps(u, t, x) = p(u(x, i:)). (In the 
general compressible case, the Eulerian PDF must be weighed on the fluctuating 
fluid density, which is tantamount to substitute pe with the Lagrangian PDF for 
the fluid parcels velocity and position 

The great advantage of the well-mixed approach, coupled with Eqn. (jH)), is 
that no knowledge of the spatio-temporal structure of the turbulent fluctuations is 
required, rather, it is the outcome, encoded in the drift coefficient, of the well mixed 
condition. This strength of the model, however, turns into weakness, whenever the 
turbulent structure plays a relevant role. A first hint that this, actually, is always 
the case, is the non-uniqueness of the solution for a given e and p_b(u) [UIB]. In the 
well mixed approach, the drift coefficient a is determined only up to a velocity curl, 
and interpretation of this freedom in terms of the properties of the flow is awkward. 

There are situations in which the structure of turbulence plays an explicit role. 
A first example is produced when coherent structures dominate the flow, a rather 
common occurrence in turbulence, which takes a dramatic form in near wall regions. 
These regions become relevant in many situations of practical interest in industrial 
flows, but also, to name a few, in the study of transport in tree canopies and in indoor 
pollution. These flows are often characterized by moderate Reynolds number, and, 
for this reason, not only the viscous scale may be not negligible, but a well developed 
inertial range may even be absent, so that the conditions justifying Eqn. Q cease 
to be valid. Now, standard techniques exist which allow for the inclusion of non- 
Gaussianity [7, and anisotropy jH E] in Lagrangian models, as well as for the effect 
of finite Reynolds numbers [H]. However, these techniques do not take into account 
the geometric structure of the turbulent fluctuations, which calls for information 
about space correlation. 

In the range of scales we are considering, another issue will come into play, 
if we are interested in modelling solid particle transport. In this case, non-tracer 
behaviours associated with inertia and gravity will begin to be felt (inertia and tra- 
jectory crossing effects jTUl HH [12] ) • This is especially true for atmospheric aerosol, 
characterized by heavy particles with relative density of the order of 1000 and par- 
ticle diameters in the range 10^^ -t- lO^/zm. Inertia effects are generally negligible in 
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ABL (atmospheric boundary layer) mesoscale modeling, but can become relevant in 
the near wall regions of wall bounded flows ^3] , which is relevant for problems of air 
conditioning and abatement of indoor pollution. Trajectory crossing effects related 
to gravity can have deep implications not only in wall bounded flows, but also for 
particulate transport in the ABL [TH] . 

These effects reflect heavily on the possibility of using Eqns. dSHHI) to model the 
Lagrangian velocity increments, and call back for the need of information on the 
spatio-temporal structure of the turbulent fluctuations. 

Some attempts in this direction were carried on in [HI [1^], but disregard of 
correlations between fluid velocity increments along solid particle trajectories lead 
to difficulties in the fluid particle limit PE] and in the implementation of the well- 
mixed condition. 

In order to be able to understand the constraints imposed by the turbulent struc- 
ture on the form of a Lagrangian model, one may try a derivation from a velocity field 
u(x, t) of prescribed statistics. If the structure function ([m*(x, t) — ^^(O, 0)]['u-' (x, t) — 
^■' (O, 0)]) scaled linearly for small space-time separations, the velocity increment be- 
tween two points lying on a trajectory would be given by: 

du= ([9t + u(x,t) ■ V]u(x,t)|u(x,t))dt + dw (4) 

with (dw) = and dw'^ = 0{dt). We could then introduce a Lagrangian model 
obeying Eqn. (^J), with 

a(v,x,t) = {[dt + u{^,t) ■ V]u(x,t)|u(x,it) = v) 
i3*^(v,x,t)dt = (dwW|u(x,t) = v) 

With the coefficients given in this equation, Eqn. would provide a Markovianized 
version of the dynamics of a particle moving in the random velocity field u(x, t). 

It turns out that, provided the structure functions for u scale linearly, the well- 
mixed technique could be imposed directly on the random field u(x, t), before any 
trajectory is defined. This means expressing the form of the conditional averages 
of the velocity derivatives: (9(U(x, t)|u(x, t)) and (Vu(x, t)|u(x, t)) in terms of the 
Eulerian PDF p£;(u, x, t). This has the important consequence that a passive tracer 
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advected by an incompressible flow, will satisfy the ergodic property 

PL(v|x,t) = ps(v,x,t) (6) 

where p2,(v|x, t) is the Lagrangian PDF for a tracer passing at the point x at time 
t, to have velocity v. Thus, the well-mixed condition imposed on the random field 
extends naturally to the Lagrangian model defined by Eqns. ((T)) and This is 
advantageous in the compressible case, where it will be shown that, contrary to 
the Thomson-87 approach, knowledge of pe is sufficient for the determination of a 
well-mixed model. 

Clearly, linear scaling at small separations does not correspond to the properties 
of a real turbulent field, which, at high Reynolds numbers is more rough, and whose 
time correlations, due to the sweep effect, have Lagrangian nature at short time 
scales This is compensated, however, by control over the large scale structure 
of the correlations, which is the relevant aspect for the determination of turbulent 
transport. 

A related issue, concerning solid particle transport, is that anomalous scaling 
of the fluid velocity increments sampled by a solid particle, are known to occur at 
sufficiently short time scales ^j. Analysis of the different ranges characterizing 
solid particle motion was carried on in JE], and Lagrangian models resolving the 
anomalous scaling range were presented in and fH], based respectively on the 
use of fractional Brownian motion and synthetic turbulence algorithms. Again, 
consideration of these short-time effects is neglected in favour of control of large 
scale geometry. 

Compared to the standard approach in Lagrangian modelling, the one proposed 
here has definite advantages. Spatio-temporal turbulent structures can be included 
in a relatively simple way. The non-uniqueness problem is solved in a simpler way, 
since only purely Eulerian properties of the flow are invoked (helicity is one example). 
The advection of passive tracers and solid particles are treated exactly on the same 
footing, hence, extension of the model to solid particle transport is automatic and 
does not need introducing additional assumptions. 

This paper is organized as follows. In section II, a local characterization of a ran- 
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dom velocity field will be given, introducing generalized "four- dimensional" Langevin 
and Fokker-Planck equations, and providing local and global existence conditions. 
The condition of local existence appears to take the form of a generalized form of 
the well-mixed condition, which will be discussed in section III. This will be used to 
calculate conditional averages in the form (Vu|u) and (9(u|u) from the property of 
dw and the Eulerian velocity PDF, and the relation with the spatio-temporal struc- 
ture of the random field will be discussed. In section IV, expressions for the noise 
amplitude (dwdw) will be derived, and their relation with the symmetric sector of 
the velocity correlation will be discussed in terms of the S0(3) technique introduced 
in jTHj. The antisymmetric sector of the velocity correlation will be discussed in sec- 
tion V, illustrating how it relates to the problem of non-uniqueness in the well-mixed 
approach, and showing how helicity and other geometrical features could be included 
in the random field. Section VI will be devoted to the derivation of a Markovian La- 
grangian model in the form of Eqn. and to presentation of its main properties. 
The relation with the Thomson-87 model [IJ will be discussed. Sections VII and 
VIII will be devoted to analysis of the Markovian approximation in the Lagrangian 
model and to proof of the ergodic property given by Eqn. Sections IX and X 
will illustrate two applications of the Lagrangian model to solid particle transport, 
respectively, in homogeneous isotropic turbulence, and in a turbulent channel flow. 
Section XI contains the conclusions. 



II. Characterization of the random velocity field 

Let us introduce a zero- mean, incompressible random velocity field u(x, t), with 
2nd order structure functions scaling linearly in the increment at small space-time 
separations. We introduce 4-vector notation: 

x^ = {x°,x^}^{t,x}, d,^^. (7) 

and stick rigorously to the Einstein convention of summation over covariant-contro- 
variant repeated indices. We have the following equation for the velocity increment: 



du' = dx^d^u' = Af,\u,x^)dx^ + dw' 



(8) 



where 

A^\u,xn = {d,u\^,t)\ui^,t)) (9) 

and (dw|u) = 0. From linearity, the contribution to the velocity structure function 
is dominated, for small values of the increments, by the correlation for dw, and we 
have: 

{dwW\u{^,t)) = {du\^,t)du\^,t)) = 0(|dx|,dt) (10) 

We limit our analysis to random velocity fields where the statistics of the velocity 
increments is independent of that of the total velocity: 

{dw'dw^\u) = (dw'dw^). (11) 

Incompressibility of the velocity field diU^ = 0, leads to the constraint, indicating 

d(Aw'Aw^) , , 

^\ = 0' = (12) 

From Eqns. (jHI) and (jlUj) . we obtain the following generalized Ito's Lemma; for a 
generic smooth function 0(u): 

d0(u) = {A^dx^ + dw')^u^^{u) + ^{dw'dw^)^u^^uJ(f){u) (13) 

and from here, we can derive an equation for the change of the 1-point PDF 
Pe{u,x'^) = p(u(x, t)), in passing from the point x^ to to the point x'^ + dx'^: 

dpE = dx'^d^pE = -du^{A;dx''pE) + ^d^^d^,{{dwW)pE) (14) 

Notice that the form of the two equations ()13|) and (fT^ is independent of incom- 
pressibility and Eqn. (fT^ . The sequence leading from Eqn. (jSj) to (fT^ is very 
suggestive, in that it generalizes the one from a Langevin to a Fokker-Planck equa- 
tion j20j. However, contrary to the case of a standard Fokker-Planck equation, 
Eqn. ()14|1 does not admit in general solution for pE- In fact, once the noise am- 
plitude (dw'^dw^) and the drift are given, Eqn. (|T4|l becomes a system of four 
partial differential equations for the single PDF pe, and this system is generally 
over- determined. In the next section, it will be shown how a generalized version of 
the well-mixed condition is able to take care of this local existence problem. 
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The ill-posedness of the problem is reflected at the global level, in the fact that 
a local definition for the "noise" increment amplitude {dw'^dw^) is not sufficient to 
define a realization for w(x, t), and consequently for u(x, t). This in contrast with 
the case of the standard Langevin equation. In fact, if we integrate Eqn. (jH}) along 
a closed curve in space-time, and consider uncorrelated increments dw along the 
curve, we will obtain in general a non-zero total velocity increment in the closed 
loop. In other words, if we disregard these correlations for w, the differential du 
entering Eqn. (jH)) will not in general be exact. 

The question becomes at this point the existence of a random velocity field with 
local structure described by Eqn. (jH)). It turns out that such a velocity field can 
be constructed explicitly, although the construction described below is by no means 
unique. 

Given a point and a direction in space time defined by the versor r'^ (r^r^ = 1), 
we can introduce the stochastic process u^{s) = 'u*(x^,r^; s) obeying the Langevin 
equation: 



The correlation functions for the stochastic process u(s), starting from the second 
order one C*-'(x^, sr^) = (u*(a;^, r^; —s/2)u^{x'^, r^; s/2)), will identify a random ve- 
locity field u{x^) whose local statistical properties are those imposed by Eqn. (jH)), 
and whose restriction to straight lines in space time will be, by construction, Marko- 
vian. As with the correlation time of the solution of a standard Langevin equation, 
the correlation length in the direction will be encoded in the drift coefficient 
r'^A^^{\i,x^). A random field realization is obtained, in the simpler Gaussian case, 
by first carrying on the principal orthogonal decomposition (POD) of C^^{x^,A^), 
and then random superposing, with the appropriate amplitudes, the resulting POD 
modes pT| . 



dff{s) = r^Afj^^u, x^)ds + dw* 

{dw'dw^)) = £{[u'{x'' + r^s) - u'{x'')][u^{xf' + r'^s) - u^ixi")]) 



s=0+ 



ds 



(15) 



8 



III. Determination of the drift 



The real meaning of Eqn. ()14|1 is to provide a consistency condition for (dw^dw^) and 
Afj^\ that could be used to generalize the Thomson-87 technique and determine from 
(dw^dw^) and pe, the expression for The difference with the standard case is 
that, instead of calculating the conditional mean (9(f*|v) of the Lagrangian velocity 
time derivative, we seek here the conditional mean of all the derivatives of the 
Eulerian velocity, namely (9^m*|u). These averages contain important information 
on the behaviour of the velocity correlation C*''(a;^, A^) = {u'^{x^ — A^/2)u^{x^ + 
A^/2)): 

C'^ix^", A^) = - A72) + R''{x^ + A'^/2)] - ^(Am'Am^') + (16) 

Here, R^^{x^) = C*-'(x^,0) indicates the Reynolds tensor, while 

C'i = hc'^ix", A^') - C'^'(x^, -A^)] (17) 

is the antisymmetric part of the velocity correlation. It is clear that the noise 
amplitude is associated with the symmetric sector of the velocity correlation, and 
for small A^: (Aw^Aw^) ~ (Am^Am^). 

Let us try to generalize the Thomson-87 approach to calculate the drift A^^ from 
Pe- It is convenient to split the drift into three pieces: 

a; = A; + — $; + — ^; (18) 

^ pE Pe 

where Aj^ is chosen to cancel the noise term in the Fokker-Planck equation p4p . 
Exploiting independence of the noise amplitude from u: 

A;dx^ = ^ {dwW)dus log PE (19) 

and Aj^, from the second of Eqn. (fT^ . is automatically traceless. The term is 
chosen to cancel the contributions to Eqn. ()14|1 from statistical non-uniformity and 
non-stationarity: 

= -d,pE (20) 
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The term \1/ is necessary to cancel the trace of $ and must be divergenceless with 
respect to u: 

dui^/ = 0, ^o' = % = (21) 

As in the Thomson-87 approach the drift is defined up to a non-unique term 
^S^* satisfying: 

El = 0, = (22) 

which, substituted into Eqn. ()14|) . will produce an identically zero contribution. 
The drift is associated with the velocity correlation through the equation: 

{u'Au') = {u'{Au^u)) = (nM/^A'^, (23) 

Let us analyze individually each of the terms in A^. Substituting Eqn. (fT^ into 
Eqn. (j23|) we see that A^ gives just the symmetric piece of the correlation, i.e. the 
— I (An* An-') in Eqn. ()23|) . The $ and \l/ terms are more easily analyzed in Fourier 
space: /(ry) = / d^ue"''''''' f{u). Using Eqns. (|TBlfW|), Eqn. ^ will read: 

{u^Au^) = -^{Aw^Aw^) - iA'^d.^i^; + ^;)|,=o (24) 

Using the fact that the generating function pE obeys pe = I — ^R^^rjirjj + 0(r]^), we 
can write, from Eqn. (j^Hj) : 

^; = ~d,R^'Vj+0{v') (25) 

so that the contribution from $ to the correlation function is, from Eqn. plj) : 
^dx^d^R^\ which accounts for the spatial inhomogeneity of the correlation [the 
^{K^^ + -R*-') term on RHS of Eqn. (fTT)|) . which is centered at ~ x^]. 

We see that A,f and account for all of the contribution to the correlations, 
which either are symmetric, or come from inhomogeneity of the statistics. We 
know at this point that both ^ ^ and S^* will be able to contribute only to the 
antisymmetric part of A^^ . We give in explicit form the contribution from 
Exploiting the first of Eqn. (pT|). we utilize the ansatz: 

i>l = i{5]7^k^^-il,^'), ^0^ = (26) 
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and, from \1>- = — $• and Eqn. (j^ : 

i^" = \{diR'') + e'^'^'r^iU + 0{rf) (27) 

with f arbitrary. The contribution to from f is traceless and can be reabsorbed 
into the non-unique term S; the final result is therefore: 

^/ = \W,{diR'^)r^, - + 0{r,') (28) 

and the contribution to the correlation function is: ^diBl'^dx'^ — diEl'^dx^], which is 
antisymmetric as required. 

Explicit expressions for the drift terms are promptly obtained in the case of 
Gaussian statistics (expressions for the case of a symmetric pE with kurtosis larger 
than three are given in the Appendix A). The velocity PDF reads therefore: 

Pe{m,x^) = Pg(u,x'^) = (87r=^||i?||)-^exp(-^5,,MV) (29) 

with Sij = {R^^)ij the Reynolds tensor inverse. In this case, the higher order terms 
in rii entering Eqns. ()25II28|) disappear and we are left with: 

A^dx" = -l{dw'dw^)Sjku'' (30) 

= \{d,K'')Suu'pE (31) 

= -^-^){diI^'')Skm + {d,R'')S,J,u^pE (32) 

We can use these explicit expressions to obtain more informations on the nature 
of the various contributions to the drift. In analogy to the case of the standard 
Langevin equation, we see that must be discontinuous at A'^ = 0. From Eqn. 
(fTTHl . discontinuity of the correlation function derivative at = is necessary to 
balance the linear scaling of d'uP' . In the coordinate system where, for the given A'^, 
is diagonal, we shall then have: 

(Am>) ~ -|A^| (33) 
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As regards we see from Eqn. ()31|) that it produces an amplification of u when 
dx^ is directed to a region in space-time where the turbulence is stronger (this is 
easy to see when oc 5"^^). 

Finally, the ^ term, turns out to produce a complicated mixture of rotations 
and amplification of the velocity vector. Indicating Ui = Siku'', and choosing the 
coordinate system so that Sji?*' = 4c(5[, and = we have from Eqn. p2|l : 

{Au^\u) = cuaA^ (AnV) = -cuiA^, {Au^\u) = 0. (34) 

If K^^ oc 6^^ , we will have Ui = Ui and the result of Eqn. (jH^ will be a rotation of u 
in the plane 12 as one moves in the direction X2- 

IV. Determination of the noise tensor 

In order to obtain the drift coefficients, which give the decay of the turbulent cor- 
relations in the various space-time directions, it is necessary first to determine the 
form of the noise tensor {dw'^dw^). In fact, it is in the noise that all the information 
on the turbulent structure is encoded (at least that part relative to the symmetric 
sector of the correlations). In the case of a Gaussian random velocity field, the noise 
tensor can be determined directly from the turbulent correlations by means of a fit 
in terms of products of exponentials with sines and cosines (a common practice in 
turbulence theory; consider e.g. the Frenkiel functions |22])- Indicating dx'^ = r^ds, 
T^r^ = 1, we fit the turbulent correlation by the expression: 

^{uHx^')uHx^' + r^s)) = CkHuHx^)u^{x^ + r'^s)), (35) 

OS 

where Cj-^ depends on the direction r'^, the mid-point position x^ + A^/2, but not 
on ds. This imposes linear dependence of the random field drift on the velocity: 

(dM>) = dx^'A; = c/mMs (36) 

(notice that Gaussian statistics, by itself, imposes linearity through the well mixed 
condition, only on the symmetric contribution to the drift A^). Using Eqn. (j36|) . 
Eqn. (j^^ takes the form: (u^du^) = Cj'^R'^ds and, from Eqn. (fT7j) . we obtain: 

{dw'dw') = {du'du^) = ^{ck'R''^ + Ck^R^') ds (37) 
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We stress that, although Eqn. (jHIj) describes the behaviors of the random field 
correlations at small separations, the coefficients c/ descend from a fit of turbulent 
correlations at finite separations. 

A very general form for the noise tensor, satisfying the incompressibility condi- 
tion d{Aw'^Aw^) I dA^ = 0, allowing association of Eqn. (jHTj) with geometric features 
of the fiow, is a superposition of terms in the form: 

{Aw'Aw^) = ^[Bi\A'') + B''J{\ - uA°))]; diB''{A) = 0. (38) 
Te 

Here, te fixes the time scales of the fiuctuations and in the Gaussian case coincides 
with the Eulerian correlation time (see next), = ^i?-, and Bl'' = lA^I^*-' + Bl'' , 
with symmetric and traceless. (For lighter notation we leave in this section the 
dependence from the space-time position unindicated) . 

We see that the presence of mixed space-time increment contributions can ac- 
count for situations in which the time correlations have Lagrangian nature. In this 
way, pure time decorrelation will take place in the reference system moving locally 
with the mean fiow u. A situation with purely Eulerian time correlation will be 
realized by putting u = 0. 

In moderately anisotropic situations, it may be expedient to expand the space 
component in spherical tensors, following the SO (3) decomposition technique 

B^\A) = J2BnA) (39) 
j=o 

where Bj indicates a combination of J-th order spherical tensors (see Appendix B). 
The symmetry of (Aw* Aw-') imposes selection rules on which spherical tensors may 
contribute; it turns out that to keep the lowest order anisotropic contribution, it is 
enough to consider spherical tensors of order J = and J = 2. The incompressibility 
condition d^iBj = gives then (the hats identify versors): 

B'^(A) = i^[(a + Ab^"'AiA^)6'^ + -(-a + (26^™ - c'")AiA„) A^A^' 
Ut 3 

-Ai[{2b^' + c'^A^' + {2b^^ + c^^)A'] + ^c"\ (40) 
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where a gives the J — part, while the tensors and , which are symmetric 
and traceless, account for the the J — 2 part. We consider next some relevant limit 
cases. 
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Isotropic turbulence 

In this case, all the spherical tensors with J > are zero. We are thus left with 
the simple expression: 

9,,2 n\A\ 1 ^ ^ 

(Aw'Aw^) = ^[|A°|5'' + -^{6'' - ttA^'A^')] (41) 

The parameter a identifies a length-scale = uxTE/a for the fiuctuations and has 
therefore the meaning of a ratio between the eddy life-time te and the eddy rotation 
time Iu/ut- 

Long axisymmetric vortices 

Let us imagine that the correlation tensor is dominated by the effect of long 
axisymmetric vortices directed along x^. Let us try to use this information to impose 
a structure to the space structure tensor B^^ defined in Eqn. ()4()|1 . Let us impose 
the condition that B'^{A) = for A = {A, 0, 0}. For A = {A, 0, 0}, we have from 
Eqn. (gOl): 

^^ = i(2a + 26" + 5cii), = a + Ab'' + Ac'', = a + Ab'' + Ac'\ 

|A| 3 ^ |A| |A| 

^ = 3c--26- ^ = 3c--26- ^ = 4c-. (42) 
|A| |A| |A| 

and we find immediately the result: 

= = 2^12 ^ 2^13 ^23 ^ Q 

3 3 

b^' = --a; b'' = b'' = -a; = --; c'' = c^^ = -. (43) 

8 ' 16 4' 8 ^ ^ 

We are free to impose the condition 6^* = c^* = for i ^ 1 and we reach the 
expression for generic A of the B components along 11 and 22: 

utB'^^A) _ 3 3X2 I 3X2 4X2 lA^A^ I ^ 

where A^ = Ag + A3 and superscripts 2 indicate squares. Analyzing Eqn. ()44|) in 
function of Ai first for A2 = and then for A2 = A_|_, it is possible to show that 
B^^ is always positive defined, as required. 
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Two-dimensional structures 

Suppose that the flow is 2- dimensional, say u = {ui,0,U3}. In this case, the 
S0(3) decomposition reduces to an S0(2) one. Keeping again only the lowest order 
anisotropic correction, we find, for A = {Ai,0, A3}: 

^^•(A) = ^[(a + 3b''^AiAm)5'' + k-a + {h'^ - c'™)AA„)AA^' 
Ut 2 

-A; [(26'^ + c'')A-'' + {2W + c'^')A^] + ?>d^] (45) 
where = 6^'^ = and the traceless condition imposes h^^ = —b^^, c^'^ = — c^^. 

Streaks 

Two-dimensional streaks along the direction of the mean fiow appear to be one of 
the characteristic structures in the viscous sublayer of wall turbulence [5J. Contrary 
to the three-dimensional case, elongated structures cannot be accommodated at 
J = 2 in an SO (2) decomposition: the resulting noise tensor would not be positive 
definite. Nonetheless, a noise expression accounting for such structures can still be 
determined. For instance, it is easy to see that, if the streaks are oriented along xi 
and the fiow is two-dimensional in the X1X3 plane, an appropriate expression for the 
noise tensor will be 

B'^{A) = a^-^SiSi (46) 
Ut 

In the non-Gaussian case, Eqn. ()36j) ceases to be valid, and the random field cor- 
relation profile ceases to be in general a simple product of exponentials and sines or 
cosines. Even if we fit the turbulent correlation with an equation like (jH^ . the ran- 
dom field correlations will not obey that equation, rather, one involving higher order 
correlations. This because of the relation, imposed by the well-mixed condition, be- 
tween non-Gaussian pE and nonlinear Af/. For instance, if we used a bi-Gaussian 
distribution to model a high kurtosis PDF [31231211, a double exponential decay 
of correlations would ensue, with the slower decay associated with the intermittent 
bursts (see the end of Appendix A) j2S]. 

For large kurtosis, the noise amplitude determines only the correlation times 
and lenghts of the fast decaying exponential. The simplest approach, in this case. 
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is to renormalize the noise amphtude, with respect to the Gaussian case, in order 
to correct for the longer correlations produced by the slowly decaying exponential. 
At the end of Appendix A, it is shown that, in order to have the desired space and 
time scales for the bursts, it is necessary to renormalize the noise amplitude by a 
factor /5 = (2/3)/c — 1 with k the kurtosis [see Eqn. (A2)]. 

V. Non-uniqueness and the antisymmetric sector 

Once the noise tensor and the PDF pE are fixed and the well-mixed condition is 
imposed, the symmetric sector of the velocity correlation is completely determined. 
The non-unique term S^* can be used to fix the structure of the anisotropic sector. 
We consider for simplicity the homogeneous case C^^{x^,A^) = C*''(A^). 

We discover immediately the following important fact: not all expressions for 
the non- unique term and consequently for A^*, lead to a statistically realizable 
C^^{Af^). This is a different face of the problem of local existence for the solutions 
of Eqn. (^H). Consider as a first example: Hj* = e^'^jU^. A contribution AC^^ = 
e^^^R^^dx^ is then added to C^^(dx), with dx = {0, 0, dx^}, that has the inadmissible 
symmetry AC*-'(dx) = — AC''*(— dx). The second example is Sq"^ = m^, Sq^ = —u^] 
in this case we find a contribution R^^dt to C^^(dt) with the inadmissible symmetry 
AC'^dt) = -AC^\dt). 

We seek a form of S^* satisfying all the required symmetries, but still sufficiently 
general to describe most geometric structures one may think of. In analogy with 
the case of the noise tensor, this can be done in the frame of an SO (3) expansion 
starting from Ca, the antisymmetric component of the correlation [see Eqn. ()17|)]. 

As with the noise, the non-unique term S^* can be determined in unique way from 
Ca in the case of Gaussian u, fitting the turbulent correlations with exponentials 
multiplying sines or cosines; in this case, S^* will depend linearly on u. Repeating 
with Ca the steps followed to obtain the noise tensor in Eqn. (jHTjl . we obtain: 

Ca = ^ [c/R''" - Cj^R'i ds (47) 
It turns out that the appropriate quantity on which to carry on the SO (3) expansion 
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is not C^, rather: 

S^kS.iC^ = ^[ci'Skj - c/5fe,]ds = r^e^eujds + ... (48) 

where Sij = {R~^)ij and the terms in the expansion indicated in the formula are 
antisymmetric spherical tensors of order J = 0,1,2 [see Eqn. (B6)]. Isolating in 
Eqn. ()36p the contribution from S^* and using Eqns. (P7j) and PHj) . we obtain 
therefore for the non-unique term: 

V = Pi^eVQ^ii?"^V + - (49) 

and it is immediate to check that the divergence free condition S^iS^* [the second 
of Eqn. (j221)] is satisfied. 

We notice that, had we carried on the expansion directly on C^, the tensor K^^ 
entering Eqn. (jl^ . and consequently the inverse Sij entering p^; in that formula 
would have been substituted by the identity matrix 6ij. The contribution E^^pe to 
the drift (here pe contains the right Sij\) would have produced therefore explosive 
behaviors (e"'"'^, a > 0) in some direction of u. What happens is that linearity of 
y4^* together with antisymmetry of C^, impose the property = cimR^^R^-' with 
Cira antisymmetric, and, keeping only the first terms in the SO (3) expansion for 
would cause loosing this property. 

We still need to enforce incompressibility, i.e. the zero trace condition S* = 
[the first of Eqn. The fact that the SO (3) expansion is carried on SikSjiC^, 

rather than on , will lead to mixing of harmonics of different order J [compare 
with the case of the noise tensor and Eqns. (B3-B4)]. It is convenient to separate 
the antisymmetric part of ^'^: 

= iL + e'rnfcC (50) 

The zero trace condition becomes: 

Si = PEii-R'^^euimu' + mi - B'lQu'} = (51) 
which must be satisfied for any uK This leads to the relation between C,^ and ^l^: 

(Ri - Rl5l)Q = ^eumm - R^f- (52) 
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The presence of the commutator [^R — -R^]^*" suggests that we should work in the 
diagonal system for K^^ . It becomes easy in this way to separate the part of which 
anticommutes with K^^ , which is simply the part out of diagonal. Solution of Eqn. 
(1^ gives then in the diagonal coordinate system, after little algebra: 

/■ _ f R22 — R33 /■ _ c R^3 — Rii /■ _ F Rn ~ R22 /p-oN 

U - ?23 p p , C,2 - Us p p , - ^2 p p • [06) 

n22 + -K33 -Kii + -rtas Kii + n22 

We are now in the position to determine the effect of the various components of Sj-' , 
on and on the velocity correlations. 

Time component 

It turns out that the /i = component of Eqn. ()49|) is associated with a combi- 
nation of rotation and strain of the velocity, as time passes, at any given position 
x; this is the Eulerian version of the mean velocity rotation along Lagrangian tra- 
jectories discussed in 0. Working in the diagonal coordinate system for K^^ , the 
contribution from will be, for instance: 

So^ = PEeoR^'u', So' = -pEeoR''^' (54) 

This turns into a pure rotation if K^^ oc 

From the point of view of SO (3), this is the J = 1 contribution to the second of 
Eqn. (B6), which is trivially symmetric in space and antisymmetric in time. 

Space component: diagonal part 

Also this component is associated with a combination of rotation and strain of 
the velocity, this time, as one moves at fixed time from one space point to another. 
Focusing e.g. on the contribution from we find in the diagonal system for R^^: 

= -pE^lR^'u^ E^^ = pE^lR^'u' (55) 

In the case R^^ oc 6^^ , this becomes a pure rotation in the plane x^x^ as one moves 
in the direction. The diagonal component of the non-unique spatial term is the 
one associated with the presence of helicity H in the turbulent field. Indicating with 
= e^-'kdjU^ the vorticity, we can write 

H = {u,uj') = {u,{uj'\u)) = e'^kiu^A/). (56) 
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Substituting the various contributions to Aj^ in the above formula, we see that the 
only terms giving non-zero result are the diagonal ones in ^. Working in the diagonal 
coordinate system, we obtain then: 



From the point of view of SO (3), this is a combination of J = 2 contributions from 
the second of Eqn. (B6) (zero trace part of ^) and J = contributions from the 
first of the same equation. 

Space component: part out of diagonal 

The effect of this component is illustrated, for the contribution from , in Figure 
^ This effect consists of a strain of the velocity components as one moves in the 



Figure 1: Sketch of the velocity lines in a coherent structure characterized by a 
non-zero value of (part out of diagonal of the non-unique spatial term). The 
velocity components in the 13 plane are arranged along strain line with expanding 
and compressing directions respectively along xi and x^. 

direction 2. We give in the equation below the non-zero matrix elements of S^* 
corresponding to ^j* (components still evaluated in the diagonal coordinate system). 



From the point of view of SO (3), this is a combination of J = 1 components from 
the third of Eqn. (B6) (the ( piece) and J = 2 components from the second of 



H = 2[^ii-R22-R33 + ^22-^11-^33 + ^33-^11-^22] 



(57) 





(58) 
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the same equation (the out of diagonal ^ piece). This is the case in which the 
incompressibihty condition needs, in order to be enforced, consideration of spherical 
tensors of different order J. 

It is possible to see that, in the non-Gaussian case, all the results obtained 
starting from Eqn. can be recovered substituting in that equation pE with 

a Gaussian PDF pc with identical S^K Notice that, if a bi-Gaussian is used to 
model pe, the ratio pc/ Pe entering the contribution to the drift will decay like a 
Gaussian for large m, when the slowly decaying Gaussian entering pE-, which decays 
slower than pc as well [see Eqns. (Al) (A2) and (A4)] become dominant. However, 
as in the case of the symmetric sector (see discussion at the end of the previous 
section), the correspondence between drift and second order correlations ceases to 
be unique as Eqn. ()36p becomes nonlinear and Eqn. ()35|) begins to involve higher 
order velocity correlations. 

VI. Derivation of Markovian Lagrangian models 

Knowing the form of the tensors A^^ and (dwMw-'), allows the derivation of La- 
grangian stochastic models. This is done most naturally setting in Eqn. (jSj) dx'^ = 
{dt, vdt}, where v is the particle velocity. This entails a Markovian assumption on 
the Lagrangian statistics, whose validity will be checked in the next two sections, 
although it is not very different from the one used in standard Lagrangian models. 

Passive tracers 

Let us write explicitly the Langevin and Fokker-Planck equations associated 
with our Lagrangian model, considering first the simpler case of a passive tracer 
dx^ = {dt, udt} where u(t) identifies the fluid velocity sampled by the moving 
particle: 

du' = ii'dt = (u ■ A* + Ao^)dt + dw' (59) 

{dt + u- V)pL + [(u ■ A^ + A^')pl] = \d^^du. (B'^pl) (60) 
where B'^ = ^^{Aw'Aw^) = ^[Bi\l) + B'^ (u)] [see Eqn. and pi(u,x,t) is 
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the PDF of finding a Lagrangian tracer at x with velocity u. 

Exactly as in the Thomson-87 approach, in the Gaussian case, the contribution 
to the drift u ■ A* + Aq^ from turbulence non-homogeneity is at most quadratic in 
u, with the quadratic terms produced by u • ($* + ^'*). However, disregarding the 
non-unique terms, the form of this contributions differs from the one discussed in 
H and p. 

Also the non-unique contribution u-H*, to lowest order in the S0(3) expansion, is 
at most quadratic in u, with the quadratic piece associated with the space component 
u ■ H*. The observation in |26j that helical contributions in Lagrangian stochastic 
models must be quadratic in the velocity is thus confirmed. 

Notice that, from the relation A = uA*^, higher orders in the S0(3) expansion 
correspond in the Lagrangian model to higher order polynomials in u contributing to 
B^^ , u ■ A* + y4g and u ■ H* + Sq*. Conversely, at the random field level, independently 
of the order in the SO (3) expansion (and for Gaussian statistics), the drift terms 
are at most linear in u. 

The important feature of the model described by Eqns. ()59ll60p is that the well- 
mixed condition imposed on the random field, encoded in Eqns. (j 1811201) and (j28|) 
[Eqns. ()3()II32|1 in the Gaussian case], translates automatically into an identical con- 
dition on the trajectories. In the incompressible case considered here, this condition 
is equivalent to the ergodic property p2.(u|x, t) = Pe{u, x, t), which will be shown to 
hold, in the next section, right thanks to the condition A] = 0. This property implies 
trivially that Eulerian averages ( ) and averages along trajectories ( )l coincide. 

At this point, the model described by Eqns. ()59II60|) . is undistinguishable from 
a model derived through the Thomson-87 technique starting from the same PDF, 
the only difference being in the form of the noise term. In its simplest form, the 
noise tensor B'^^ is isotropic, and is obtained by setting A = uA° in Eqn. ()41|) and 
deriving with respect to A": 

= + ^(^s'J _ i^)]. (61) 

te Ut 3|u|^ 

This expression must be compared with the one in the Thomson-87 approach: B'^^ = 
(5'-'Coe. (As a technical aside, notice that we started in section II with an additive 
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noise and we have arrived here at a multiphcative noise term, which is automatically 
intended, in the approach that we have followed, in the Ito sense |2Uj). 

The difference in the analytical expressions underlies a difference in physical 
interpretation: while in the Thomson-87 technique, B^^dt is precisely the Lagrangian 
time structure function for inertial time separation, in our approach, it is a non- 
universal quantity whose form is determined in function of the large scale turbulence 
geometry. In the Thomson-87 approach, the time scale is fixed by the viscous 
dissipation e, which fixes the expression for the Lagrangian correlation time 

In our approach, the time scale is fixed directly by te- To be precise, the association 
between r/, and e in the Thomson-87 model, is strictly valid only in the Gaussian 
case |2^j. Also in our approach, however, non-Gaussian statistics leads to te not 
being directly associated with the Eulerian time scales, but only with the fast part 
of the correlation decay (see end of Appendix A). 

The two approaches depend on dimensionless constants, Co and a, which can 
be related in semi-quantitative way. As discussed in correspondence to Eqn. (jH)), 
the parameter a identifies a characteristic length for the random field lu = utTe/o,, 
which, at least in the Gaussian case, corresponds to the integral lenght of the turbu- 
lence. Substiting the estimate for the viscous dissipation from our model e ~ Uj^/lu 
into Eqn. and setting Te ~ t^, we obtain: 

a ~ Cq^ 

This tells us that the Thomson-87 model cannot be recovered from Eqn. for 



2u- — 

finite Cq, by setting simultaneously a = and — - = Cot. The equivalent limits 
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TE 



a ^ and Co ^ oo correspond to the regime of te much shorter than the eddy 
rotation time corresponding to the Kraichnan model j2Hl- The way in which this 
limit is carried out, however, is different in the two approaches: in ours, it is the 
turbulence integral scale that is sent to infinity [see comment after Eqn. (PT|) ]: in 
the standard approach, the diffusive limit Cq oo corresponds directly to the one 

TE-^0. 
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Solid particles 

Let us pass to the analysis of the sohd particle case. The solid particle dynamics 
obeys an equation, which in quite general form (neglecting memory effects associated 
with the Basset force) can be written as: 

V = F(v,u) (63) 

A general form was derived by [HD] , who neglect lift effects • These equations are 
all derived in the limit of particle diameter d small with respect to the scales of the 
flow, and low particle Reynolds number Rcp = d\u — v|/z/, where u is the kinematic 
viscosity. We will consider Eqn. ()63|) in simplified form by accounting only for linear 
Stokes drag and gravity: 

V = "-" + "^ (64) 

where ts = 1/18 'yd'^/u, with 7 the density of the particle relative to that of the 
fluid, is the Stokes time, and is the particle terminal velocity in a uniform force 
field and a quiescent fluid. In the case of gravity: = Tsg with g the gravitational 
acceleration. More in general, vg, may account for body forces like the effect of the 
Saffman lift [31j. 

The analogue of Eqn. (j^Hl), in the solid particle case will be obtained putting in 
Eqn. (jHl) dx^ = {dt, vdt}, with v the solid particle velocity: 

du' = u'dt = (v ■ + A^;)dt + dw' (65) 

where now u(t) is the fluid velocity sampled by the solid particle and {dw^dw^) = 
B'^dt = ^[Bi\l) + B'^{-v)]dt. Notice that the drift tensor A^' still depends on u, 
while (dw^dw^) depends only on v. The Lagrangian PDF pl(u, v, x, t) will obey the 
Fokker-Planck equation: 

(5i + V • V)pi + d,. (F'pl) + d^. [(v ■ A'^ + AoOpl] = \d.^d^, (B'^pl) (66) 

Equations (|63|) . (j65j) and (j66j) are in the standard form for a "two-fluid" Lagrangian 
model for solid particle transport, i.e. a model in which the fluid and solid phase 
are taken into account at the same time and are treated on the same footing. All 
problems in the fluid limit, present in models in which the separation of fluid and 
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solid particle trajectories was considered without accounting for the geometry of the 
process ^H] are clearly avoided: when inertia and gravity are sent to zero, the fluid 
case described by Eqns. (I59II60|1 is automatically recovered. 

We can easily estimate the turbophoretic drift, i.e. the component of particle 
transport due to the turbulence intensity gradient [321 EBj- Multiplying Eqn. by 
f * and integrating in v and u we obtain at stationary state and for uniform concen- 
tration dj{v^v^)i^ = (F*(v,u))l, which, for linear F*, can be inverted to obtain (f*)L 
[subscript L indicates that we are averaging over the Lagrangian PDF Pl(u, v, x, t)]. 
In the Stokesian case described by Eqn. (j66|) . and small Stokes number St = ts/te, 
we can approximate {v^v^)l = R^-' , and we obtain: 

{v\ = -Tsd.R'' + vh (67) 

We can try to understand Eqns. ()65ll67p from the point of view of a model satisfying 
the well mixed condition. The particle flow, due to the effect of inertia, is compress- 
ible and preferential concentration phenomena are known to occur [SH]. Therefore, 
we do not expect in general the ergodic property pi(u, v|x, t) = p£;(v,u, x, t) to 
be satisfied. Turbophoresis provides the simplest illustration of this phenomenon. 
Averaging Eqn. (j64|) over pi(u, v,x, t) and combining with Eqn. (j67|) . we obtain in 
fact the relation: 

{u\ = -rsd.R^ ^ {u') = 

i.e., in inhomogeneous turbulence conditions, Eulerian and Lagrangian averages give 
different results. 

For these reasons, in the Thomson-87 approach, a two-fluid solid particle trans- 
port model, would require knowledge, in some reference situation, of the Lagrangian 
PDF pl(u, v,x, t), meaning that informations must be available on both the mean 
particle concentration 6'(x, t) = J d^Md^fpL(u, v, x, t) and the conditional PDF (the 
PDF along a single particle trajectory) p2,(u, v|x, t). Notice that this may imply 
substituting Eqn. or with a model equation whose coefficient are deter- 
mined by the well mixed condition. Actually, one- fluid models exist, in which only 
the particle phase is considered and only the PDF pi(v|x, t) has to be known, while 
6 is obtained from the continuity equation dfO + V ■ ((v)l^) = 39 . 
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In our approach, knowledge of pe is sufficient to determine the form of the 
equation for u (j65j) (the one for v is unchanged) without any assumption on the 
form of Pl(u, v, x, t). This will be shown in Section IX to produce automatically the 
correct form of the Lagrangian correlations for the fluid velocity along solid particle 
trajectories, accounting for the effect of inertia and trajectory crossing [TT]. Notice 
that, imposing the well mixed condition on pl(u|x, t) = j d^fpL(u, v|x, t), within 
a simplifying ergodic hypothesis pi = Pe, is not sufficient to obtain these correct 
behaviors; the anisotropic renormalization of the correlation times may be accounted 
for only by ad-hoc modification of the expression for the noise tensor Eqn. © p]. 

VII. Lagrangian one-time statistics and ergodic 
properties 

As discussed at the end of Section II, we can imagine Eqn. (jHl) as giving the local 
behaviour of a random velocity field whose restriction to straight lines in space-time 
are Markovian processes. This allowed to have a random velocity field with Eulerian 
correlations in time and space, which are both well defined and easy to calculate. 
Unfortunately, unless the a limit of the Kraichnan model is considered [2H], 
it is not possible to hypothesize at the same time a Markovian behaviour along 
trajectories. In consequence of this, the Lagrangian statistics becomes a complicated 
business. 

However, it turns out that different statistical quantities are affected by the 
presence of memory in qualitative different ways and there are situations in which 
Markovianization of the trajectories becomes appropriate. Let us try to understand 
what happens in detail. 

The central quantity one needs for a description of Lagrangian statistics are 
conditional probabilities in the form 



where = X(tfe) = {u(x(tfc), tfc), x(tfc)}, k = 0,...n. Let us consider for now the 
simplest case of a passive tracer. Such conditional probabilities could be obtained 
ideally by carrying on a Montecarlo of trajectories originating from {t„,X„} and 




(68) 
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sampling the particle positions and velocities at t = tk, k = n — 1, ...1. Let us focus 
on the case n = 1, which presents already all the difficulties due to memory. Actually, 
the transition probability p(Xo|Xi) gives precisely the evolution of a cloud of tracers 
from an instantaneous release, i.e. a puff; from the same transition probability, also 
the Lagrangian correlation time could be determined and the calculation will be 
illustrated in the next section. Suppose we have a set of trajectories starting at time 
ti with initial condition Xi, whose form is known up to time t. This allows us to 
reconstruct pL(X(t)|Xi). The conditional probability at the instant t + dt will be 
given by the formula: 

p^(X(t + dt)|Xi) 

= J d6X(t)pL(X(t + dt)|X(t),Xi)pi(X(t)|Xi) (69) 

which corresponds to first summing all the trajectories going from {ti,^i} to {t + 
dt, X(t + dt)} passing through {t,X(t)}, and then summing over X(t). Now, to 
determine p(X(t+dt)|X(t), Xi), we could average first on the part of the trajectories 
going from {ti,Xi} to {t,X(t)} and then on that going from {t, X(t)} to {t + 
dt, X(t + dt)}. From the point of view of a Montecarlo, this means that we can 
consider an ensemble of fictitious trajectories whose dynamics is only conditioned 
to the initial condition {ti,Xi} and to the current position {t,X(t)}. 

We thus reach the not so obvious conclusion that, to determine the evolution of 
a PDF with conditions at n previous instants, we need to study a dynamics condi- 
tioned to n + 1 instants, but we do not need the whole trajectory history. Because 
of this, if we are interested in a 1-time PDF, Markovianization of the dynamics will 
be an appropriate procedure. 

This fact allows us to verify analytically that the 1-point velocity PDF sampled 
by a passive tracer coincides with the Eulerian PDF; in other words the ergodic 
property one expects from incompressibility is satisfied. 

In the case of a passive tracer, dx'^ = {dt, udt} where u(t) identifies the fluid ve- 
locity sampled by the moving particle, and Eqns. and (jUUj) will be the Langevin 
and Fokker-Planck equations associated with the Markovianized dynamics. From 
incompressibility [see Eqn. fjl2|) ]. and from the properties of the drift components 
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A, $ and \1/ [see Eqns. P8II21|) ]. we can write: 

{dt + n- V)pl + + = -d^. [{A/u' + A,' - ]^B'^du,)pL] (70) 

Setting Pl = Pe, from Eqns. (fT^ and the time component of Eqn. (pUj). Eqn. (fTOj) 
reduces to u'diPE = —u^du^^/, which, from Eqn. ((20)), is an identity. The ergodic 
property is thus satisfied. 

This is not a trivial property. We can easily construct a counter-example in 
which the incompressibility property d/\t{Aw'^Aw^) = [the second of Eqn. (fT2|) ] 
is not satisfied and ergodicity is violated. Considering for simplicity stationary 
homogeneous turbulence (hence <I>^' = ^^^^ = 0) and choosing S^* = 0, we have in 
fact, setting in Eqn. (fT^ dx'' = {dt,udt}: 

n-A' + Ao' = ^B'W^AogpE (71) 

while Eqn. (jUnj) dictates: 

u-A^ + A,' = ^B'^d^, hgpL + \d^.&^ (72) 

Combining Eqns. (f?^ and (fTSj) . leads to a differential equation for pL which, due 
to homogeneity of i3 in |u|, can be integrated along the direction u = u/|u|: 

Pl(u) = const pEivi) exp \ ~ I ds Ub'^) jid.^kB^^] \ (73) 

K J \- J u=us J 

Taking a noise term not satisfying incompressibility, e.g. i?*-'(A) = ^^lAj^*-', we 
would obtain Pl(u) = const |u|^^p£;(u) and ergodicity violation. 

We can repeat the calculation to check for departures from ergodicity in the 
solid particle case. Ergodicity means in this case that the fluid velocity distribution 
sampled by the solid particle 
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'u|x) = j p/^(u, v|x)(iv (74) 



coincides with p£;(u, x). In all the Montecarlo simulations that we have carried on, 
described in detail in Section IX, we have found that, despite compressibility of the 
solid particle flow, ergodicity was satisfied in isotropic homogeneous conditions. The 
mechanism seems to be the following. 
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In homogeneous stationary conditions, the Fokker-Planck equation for the dis- 
tribution pi(u, v) will read, from Eqn. (jUUI) : 

d^F'pL) + 9n4(v ■ + v4o>l] = ^d^^d^^B'^pL (75) 

where B^^ = ;B*-'(v). Exploiting well- mixed [see Eqn. (jl9j) ] and setting from isotropy 
2 = 0, this equation can be rewritten in the form: 

dAF'pL) + d^.[^B'^{d^,pL - PlOu. \ogpE)] = (76) 

and, integrating in d^v, we reach the following equation for the deviation from 
ergodicity: 

9„,[pi((i3*>)9„. hgpJpE + du,{B'^u))] = (77) 

We see that non-ergodic behaviours are associated with the divergence of the average 
over solid particle trajectories of the velocity structure function: 

du.{B'^n) = j d\9„.pL(v|u)S^^(v) (78) 

In the case of solid particles, for which {B^^\\i) ^ i3*''(u), we would expect in general 
|u) 7^ 0. This turns out not to be true, however, when turbulence is isotropic. 
Let us show how this happens. 

We can decompose VuPl(v|u) in spherical vectors depending on v [see Eqns. 
(B7-B8)]: 

VnPL(v|u) = poiv + pii(u ■ v)v + P12U + ... (79) 

where, from isotropy, pik = pifc(|v|, |u|). Higher harmonics (not indicated) are by 
construction orthogonal (see Appendix B). From Eqn. (PT|) . we see that only the 
term 

h = pii(u ■ v)u + P12U (80) 

can contribute to d^j (B^^u). In order for this contribution to be zero, it is sufficient 
that the curl with respect to v of h be identically zero: 

[|v|"^(9|v|Pl2 - Pii]v X u = (81) 
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so that h can be written in the form of a potential term V^(7(u, v), and, substituting 
into Eqn. (f7S|) and integrating by parts: 

= -j d\(7(u,v)9,.i3^^ (v) = 0. 

Now, we can obtain Eqn. (jHTjl simply by imposing the condition, from isotropy: 

a„,(|v|W|u) = 9„.(|vrt;V|u), (82) 

with i ^ j ^ k and n > 0. In fact, writing the averages in explicit form, Eqn. (jH^ 
can be shown to be equivalent to: 

J d=^t;|vr+V[9,.9,, - 9,.9„.]Pl(v|u) (83) 

and again, from orthogonality of the decomposition, only the pu and pi2 terms in 
VmPl could contribute. Hence, exploiting the fact that pij depends only on |u| and 
|v|, Eqn. (jHHj) becomes equivalent to 

j d3t;[|v|-i9HPi2-pii]|vr+i = 

which implies Eqn. (|HTj) and satisfaction of the ergodic property. 

VIII. Two-time statistics and the Lagrangian cor- 
relation time 

Explicit determination of the Lagrangian dynamics taking into account memory of 
an initial condition is possible when the u(x, t) is isotropic, homogeneous and Gaus- 
sian. It thus becomes possible to estimate the error implied in the Markovianization 
of the trajectories. The simplest estimator is the Lagrangian correlation time 

1 f°° 

= ^ I dt(u(t) ■ u(0)) (84) 

where u{t) = u(x(t),t) and we are considering passive tracers. As discussed at 
the start of the previous section, we need an evolution equation for the trajectory 
{u(t),x(t)}, given an initial condition at t = [for simplicity, fix x(0) = 0]. The 
starting point is the following decomposition for the tracer velocity: 

u(t + A) = (u(t + A)|u(t),x(t);u(0),0) + Aw (85) 
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plus knowledge of the conditional averages: 

(u(t + A)|u(t),x(t);u(0),0) and (AwAw|u(t), x(t); u(0), 0) (86) 

As discussed in detail in Appendix C, these averages can be obtained from the 
correlation between velocities at {0,0}, {t,x(t)} and {t + A, x()f: + A)}. We identify 
correlations between points on a trajectory by: 

C''J{^) = {u\t)ui{t + ^)), (87) 

If u is Gaussian, homogeneous and isotropic, these correlations can be expressed in 
analytical form. The mean rate of fluid velocity change along a generic space-time 
direction {1, V} will in this case take the form: 



te 



|V|\ . / 2|V| 
1 + a—]u\ + 1 + ]u\ 



(88) 



where and uy are the components of u perpendicular and parallel to the fixed 
direction V, and we have used Eqns. (j^ . (IHUj) and (HH), and the expression R^^ = 

Solving the equations for the correlation function along {1, V}: 

^(«^(Vt,t)M^(0,0)) = {[A,' + F^A/]m^(0,0)) (89) 
and introducing longitudinal and transverse projectors 

ny(v) = ^, nY(v) = 5^^-n;f(v) (90) 



we obtain: 



where 



(w^(vt,t)w^(o,o)) = n^i(v)Cy^(t) + nj(v)Cy||(t) (91) 

Cy±(t) =M|expf- fl + ^)— ) (92) 
V V Ut ' TrJ 



and 



Cv,(i) = 4-p(-(l + a|^);^) (93) 

We shall need also the inverse Di^iti — tm) of the correlation matrix {u''(ti)u^ (tm)), 
l,m = 1,2] ti = 0, t2 = t, defined by the relation: 

Y^D^.iU - t„,){u'{t^)u\K)) = Y,{^\Qu'itrn))D,,{t„, - ti) = 6^6in (94) 
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From Eqns. ()91II93|) . we find: 

Dijiti - tm) = nJ(U)Dx(tz - tm) + n;i.(U)D||(t, - (95) 

where U = u(t) — u(0), 

D 1 ( Cu±iO) -Cu±{t)\ 

c'^^o) - ci,^{t) y-CuAt) Cum) ^ ' 

and we have similar expression for D\\. At this point, we can obtain from Eqn. 
(C7) the expression for the average evolution of the velocity along a trajectory, 
conditioned to an initial condition at time zero: 

{u\t + A)|u(t), x(t); u(0), 0) = [&^{A)Djk{t) + C'^{t + A)D,k{0)h\0) 



+ [C*^(A)D,fe(0) + C'^{t + A)D^k{t)V{t) (97) 

As obvious, memory of the initial condition at time zero is lost when t ^ oo. Notice 
that, if all points {0, 0}, {t, x(t)} and {t + A, x(t + A)} are aligned along the same 
space-time direction {1, V}, all the and entering the correlation functions in 
Eqn. ()95|) will project along or perpendicular the same vector V. In this case the 
components of u parallel and perpendicular to V decouple and the indices disappear; 
for instance: 

{u\\{t + A)|u(t),x(t); u(0), 0) = [Cv\\{A)Dv\\{t) + Cv\\{t + A)Dy||(0)]n||(0) 



+ [Cv\\ {A)Dv\\ (0) + Cv\\ {t + A)Dv\\ {t)]u\\ (t) (98) 

and all the correlations have the same decay rate fixed by Eqn. (jU^ . It is then easy 
to show that the first term on the RHS of Eqn. (j98|) disappears and we have 

(M||(t + A)|u(t),x(t);u(0),0) =M||(t)expf- fl + a^)— ) (99) 

Hence {u\\{t + A)|u(t), x(t); u(0), 0) = {u\\{t + A)|u(t), x(t)). If the trajectory is 
developing along a straight line, we will recover Markovian statistics as required. 
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Expressing Aw as the difference between {u(t + A)|u(t), x(t); u(0), 0) and u(t + 
A), and substituting into Eqn. (C8), we obtain instead, for the ffuctuation term: 

{Aw'Aw^\u{t),x{t);u{0),0) = 

= C'^iO) - Dki{0)[C^\A)C''\A) + &'{t + A)C''^{t + A)] 



-Dkiit)[C^\A)C''^{t + A) + C^'{t + A)C'=^'(A)] 



(100) 



In Figure El we compare the result of a Montecarlo for the Lagrangian correlation 
time Eqn. (jH^ using the exact dynamics described by Eqns. ^7\i and pOOp . 

with that obtained from the Markovianized version given by Eqn. (jSHI)- As could 
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Figure 2: Dependence of the Lagrangian correlation time on the ratio a between 
eddy life time and eddy rotation time; □ exact; O Markovian approximation. 



be guessed, the Markovian approximation becomes exact in the a = limit, when 
the trajectory, in a correlation time, remains close to the time line {1,0}. At least 
in this case, the choice given by Eqn. (jisp . of Markovian statistics along rectilinear 
cuts in space-time, is the most appropriate. 
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IX. Solid particle transport in homogeneous isotropic 
turbulence 



Inertia and crossing trajectory effects determine a substantial change in the statistics 
of fluid velocities sampled by the solid particle with respect to that of passive tracer 
velocities. 

Several authors reserved particular attention to the long time behaviour of cor- 
relation functions of sampled fluid velocities and long-time particle diffusion coeffi- 
cients unmiiEiiEni- 

Following Csanady jTU], Sawford and Guest [12] kept into account the effect of 
gravity produced trajectory crossing by a suitable assumption on the renormalization 
of the correlation time of the fluid velocity sampled by the falling particle. Their 
model was applied to grid-generated turbulence and their results were found to 
agree with experimental wind tunnel data. It is not clear, however, how much this 
approach can be extended to generic non-homogeneous and non-stationary turbulent 
flows, especially in the case of strong turbulence gradients jSEI- 

Free-flight models (see also for a brief review) are known to make un- 
physical assumptions about the velocity that particles assume when they are pro- 
jected towards the wall from the buffer and logarithmic regions. As regards the 
eddy-interaction model of Kallio and Reeks , this has been shown in [38J not to 
satisfy the well-mixed condition. The model described in [38J improved this aspect, 
but without reproducing the build-up of concentration. The issue, to be discussed 
in the next section, is the difficulty in isolating near wall solid particle accumula- 
tion effects from spurious concentration build-up from unproper treatment of the 
well-mixed condition. 

A recent advance was obtained in jSH], but in this approach a turbophoretic force 
had to be introduced from the outside, whereas in our approach the turbophoretic 
flux results in self-consistent way from the dynamics [see Eqn. ()67|)]. 

The central role in the solid particle dispersion is played by the correlation time 
t'j^, of the fluid velocity sampled by the solid particles. In particular, with r^(||) 
we indicate the longitudinal effective Lagrangian time, i.e., along the direction of 
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gravity, and with Tj^(-L) the transverse effective Lagrangian time. 

Let us brieffy summarize the main properties of T^(i). When gravity is dominant 
{vg ^ ut), the correlation function of sampled fluid velocities decays faster than 
that of passive tracers and t^(«) < ti, where i is referred to longitudinal (||) or 
transverse (±) [121111]. Taking vq = QTs with g fixed, we see that r'j^ii) decreases 
from the value r^, corresponding to ts = 0, to zero as ts increases. Furthermore, 
the correlation functions do not decay in the same way in all directions. Due to the 
continuity effect described in ^0], the decay is slower in the direction of gravity, so 
that the longitudinal Lagrangian time t^(||) is longer than the transverse one r^(_L). 

In the inertia dominated case {vq <^ mt) , the sampled correlation function decays 
slower than that of passive tracer velocities and > tl. The limit r5 ^ is the 
same as above, but now rj^ increases with ts and, in the limit ts — * oo, tends to the 
Eulerian time scale te jTH I34j . 

We will show shortly how all these effects are automatically reproduced in our 
approach. 

We consider a Gaussian homogeneous and stationary isotropic zero-mean random 
velocity field. Thus, the drift term is given by Eqn. ()3()|1 with S'*-' = 5^^ /u\, the 
PDF is given by Eqn. (j^^ and the noise tensor is isotropic [see Eqn. Mlj) ]. The 
terms and ^ ^ are zero for homogeneity and the non-unique terms S^* are zero 
for isotropy. 

From now on in this section, we rewrite the equations expressing the velocities 
u and V and time t in units of ut and te respectively. Note that, in this way, 
Tg becomes equivalent to the Stokes number St jSU], i-e., the ratio between the 
Stokes time and a flow time scale {je in this case). As regards gravity, it results 
vg = St/Fr, being Fr = g te/ut the Froude number related to the magnitude of 
the gravity g with respect to turbulence scales 

Under these conditions, Eqn. ()65|) for the sampled fluid velocity u{t) and the 
associated expression for the noise tensor B''^ will take the simplified form: 

jdu' = + a \v\)u^dt + ^v^UjV'dt + dw' 

[{dw'dw^) =2[{l + a \^r\)6'^ - !|\^r\ v'v^]dt ' 
with V = v/|v| and v the particle velocity, whose dynamics is given by Eqn. (jM|) . 
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Applying to Eqn. ()101|) the projectors defined in Eqn. ()90|). with V = v, it is easily 
seen that du is split into a longitudinal and a transverse component, characterized 
by different values of the drift: 1 + 2/3 a|v| for the fluid velocity component parallel 
to the particle velocity and 1 + a|v| for the normal component. As the role of 
gravity increases, the symmetry breaking of particle motion due to the presence of a 
preferential direction, i.e., the direction of gravity, involves a separation between the 
longitudinal and transverse time scale (continuity effect). In the gravity dominated 
case, vg ^ 1, we have |v| ~ vg, with the result: 

' r'M^^^. (102) 



^^"^ 2avG' ' avG ri(±) 2 

It is difficult to obtain an analytical solution for the PDF pi(u, v) and the velocity 
correlation functions {u^{t)u^{0)) because of the multiplicative noise. For this reason, 
numerical simulations by means of Montecarlo technique have been performed to 
obtain solutions of Eqns. (jMjl and (jlOlll . Following Yeung and Pope [101 (see also 
|41j). we choose a value t^/te ~ 0.5, which, from FigureEl corresponds to a = 0.65. 

As mentioned in Section VII, the ergodic property has been verified to hold also 
in the solid particle case. Both in the case of Gaussian statistics and of an isotropic 
kurtosis described by Eqn. (A3) (see Appendix A), with and without gravity, the 
marginal Lagrangian PDF Pl(u) defined in Eqn. ()74|1 has been found to coincide, 
to within numerical error, with the Eulerian PDF Pe{u). 

In the presence of gravity, this means that the average of the sampled fluid 
velocity (u)^, coincides with its Eulerian counterpart (u), which is zero, and that, 
therefore, no renormalization is produced on the value of the terminal velocity Vfj. 

Another non trivial result from the numerical simulation is that the correlation 
functions of both passive tracers velocities {ts = and vg = 0) and of sampled 
fluid velocities appears to decay exponentially as in the Gaussian case for standard 
Lagrangian models. 

In Figure El the effective transverse Lagrangian times Tj^(-L) have been plotted 
as function of Ts ( in units of Te) for different values of Vg (the Lagrangian time 
scale of passive tracer has been reported for a comparison). In the range fc < 1 
(inertia dominant) the curves are increasing from tl/te = 0.52 and tend to about 
1 as r5 ^ oo (i.e., the Lagrangian time tends to the Eulerian time). For f ^ > 1 
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Figure 3: Behaviour of r^(-L) as a function of the Stokes time ts for different values of 
the adimensional terminal velocity vq (dimensionless units). O vg = 0; + vg = 0.1; 
n vg = 0.5; X fG = 1; A = 2; * = 5. Reference line at = t/, = 0.52. 

(gravity dominant) the curves loose their dependence on Ts and the correlation time 
is approximately equal to the eddy crossing time given (always in dimensionless 
units) by Vq^. This is in agreement with the asymptotic formulae in Eqns. (jlU2j) . 
Figure m shows the behaviour of t^(||) and t{^(-L) as function of vg for fixed r^. As 
expected, each curve tends to a constant value in the limit vg ^ 0. The longitudinal 
times are always longer than the transverse ones and, in the limit fG — ^ oo, they 
collapse onto two different curves, whose ratio is 3/2 as predicted by Eqn. (jl(J2j) . 

An exponential decay of the sampled fluid velocity correlation function allows 
easy analytical calculation of the correlation function for v jSH |H3 The last 
one reads, for i = ||, ±: 



■t/rs 



(103) 



+ 



where 




(104) 
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Figure 4: Behaviour of r^(||) and t];^(-L) as a function of the terminal velocity vg = 
for different values of ts (dimensionless units). Longitudinal: * = 10; x ts = 1] 
+ Ts = 0.1. Transverse: A ts = 10; □ = 1; O = 0.1. 

By using Taylor's theorem, the (long-time) diffusion coefficients 

= J lim 7((a;,(t) - Xi{0)f), i= ||, _L 

can be expressed in terms of the Lagrangian correlation times for v as follows 

t » r,(z) : = r,(z)C;(0) = T'^{t) 

Thus, the diffusion coefficients K,{i) will behave exactly as the effective Lagrangian 
times T^(i). The adimensional diffusion coefficient of passive tracers is simply given 
hj K, = tl = 0.52. Hence, in agreement with ^1], K{i) will be larger than in the 
passive scalar case when inertia is dominant, smaller when gravity is dominant. 
Furthermore, when gravity is dominant, the longitudinal solid particle diffusion 
coefficient will be larger than the transverse one. 

X. Solid particle transport in turbulent channel 
flow 

We focus in this section on phenomena of accumulation and deposition associated 
with the interaction of inertial particles with the inhomogeneity of the flow and 
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the presence of solid boundaries. Starting from the work of McLaughhn the 
reference situation that is typically considered, both to identify the main features of 
particle transport, and to test the functionality of transport models, is that of the 
turbulent channel flow. 

We have tested our model in its simplest form, with a Gaussian PDF, an isotropic 
noise and non-unique term S^* set to zero. We recall that, in this form, the model 
is described by Eqns. fj64ll65p . with the drift given by Eqns. (fTHj) and ()30ll32p . and 
the noise by Eqn. (jlT|) with A = vA° = vAt. As in the homogeneous-isotropic 
turbulence case, we have set the free parameter a = 0.65, corresponding to the value 
of the ratio between Lagrangian and Eulerian correlation time t^/te = 0.52 j40|, 1^ . 
We adopt standard wall variables identified where necessary with +, normalized 
with the friction velocity and the reference length and time scales X2 = z//m* and 
r^, = Xg/w* where u is the kinematic viscosity. 

For the Lagrangian correlation time we have used the interpolation formula: 



where X2 identifies the cross-stream direction (we take xi and X3 respectively in 
the streamwise and spanwise direction). For < 140, Eqn. ()105p coincides with 
the interpolation formula quoted in [3T. Thus, Eulerian time scales span from 
^max ^ ]^^g ^^^q channel centre to r™" ~ 14 at the walls. 

We have considered neither the effect of gravity, nor that of Brownian motion. 
The second may be important in the case of sub-micrometer particles. We have 
included, instead, the contribution from the Saffman lift indicating as usual 
with d the particle diameter and 7 the particle to fluid density ratio: 



which is known to contribute to the solid particle dynamics in the range < 10 



As input data for the model, we have utilized the 1-point statistics from the DNS 
by Kim, Moin and Moser The channel width Lc is nearly 360 wall units and 



tl = 7.122 + 0.5731 xt - 0.00129(x^)2 
= -19.902 + .959 x^ - .00267(x^)2 



xt < 140 

140 <xt < 180 



(105) 




du 



du 
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the Reynolds number is of order 3300, based on the maximum mean velocity and 
the channel half-width. 

We have carried on Montecarlo simulations with Ntot = 10000 particles, uni- 
formly distributed at the initial time. In order to obtain, after a suitable time, 
a stationary concentration profile, we have introduced a source of particles at the 
channel centre to balance the deposition flux at the walls. As in for these con- 
ditions and for quite all Stokes times, it seems that a simulation time T^jm of about 
700 is sufficient to achieve a stationary distribution for the solid particles. 

As a validation. Figure El shows that the model reproduces, in the fluid particle 
case, the input statistics. Furthermore, the well-mixed condition has been verified: 
despite the possible numerical complications arising from the presence of a multi- 
plicative noise, a uniform passive tracer concentration profile is preserved in time 
and no tracer deposition on the walls takes place. In fig. IHl we give the profile of 

3 
2.5 

2 
1.5 

R'^ 1 

0.5 



-0.5 

-1 

50 100 150 200 250 300 350 
^2 

Figure 5: Comparison between input statistics for the Reynolds tensor K^^ , from 
DNS jlH] , and simulated data from Montecarlo. The data are almost undistinguish- 
able: (a) (b) R^^, (c) R^^, (d) R^'\ 

the fiuctuation amplitude for the normal velocity of a particle with ts = 60. The 
Montecarlo data strongly differ from the profile obtained from the homogeneous 
isotropic turbulence estimate provided by Eqn. (jl04p and illustrate the difficulty 
in the a-priori determination of a reference PDF pi(v,x, t) in one-fiuid Lagrangian 
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Figure 6: Comparison between the Montecarlo simulation results for {v"^) profile 
(a), and the homogeneous turbulence estimate for the same quantity (b). 

models (see discussion at the end of Section VI). 

In Fig. [7| we give account of the particle concentration build-up in the near 
wall regions. The peak height appears to increase with ts up to r5 ^ 10 and to 
decrease afterwards; the same decrease was observed in jlHI- In agreement with 
both 03] and jlH], and in contrast with the one- fluid model in (SHj, we observe that 
the concentration maximum occurs in the viscous sublayer at < 1. Conversely, 
numerical data on the peak height present in literature show a definite scatter; 
anyway, our data are closer to those in [33] than in [31], with an over-estimation of 
the order of 50% with respect to the first. 

As regards particle deposition, we have studied the dependence on ts of the 
deposition flux 

~ T N 

-'- Sim-'- ^tot 

being the channel width, A'^^ the number of deposited particles in the simulation 
time Tsim and Ntot the total number of particle simultaneously present in the channel. 
In our simulations, we consider a particle deposited, when its distance from a wall 
is smaller than its radius d/2. Assuming that air is the suspending medium, we fix 
for the density ratio the value 7 = 1000, and for the viscosity v = 0.15cm^/s; from 
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Figure 7: Concentration profile 6 vs. x^. (a) ts = 3; (b) ts = 10; (c) ts = 20; (d) 
Ts = 30; (e) ts = 60. 

relation ts = 1/18'yd'^ /u, the particle diameter will then be, in wall units: 

d ~ 0.134r| 

Our results are shown in Fig. |H1 and compared with experimental data by jTTl |3H1 , 
and with an example of one- fluid Lagrangian model ^\ . The agreement is good with 
respect to the data in [48 , apart of a slight over-estimation in the range ts > 10. 
On the contrary, our model performs much better than the one-fluid model in the 
range < 10. As expected, the effect of the Saffman lift is felt only in the range 
Ts < 10; in any case, the contribution to both deposition and particle accumulation 
appears to be small. 

XI. Conclusions 

We have studied the statistical properties of trajectories extracted from a random 
velocity field with non-zero correlation time, analyzing the conditions for a Marko- 
vian approximation of the Lagrangian velocity. Our result is that a generalized form 
of the Thomson-87 well-mixed condition [3] can be derived also in the case of random 
fields, provided their velocity structure functions scale linearly at small space-time 
separations. In the incompressible case, the Markovian approximation for the La- 
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Figure 8: Comparison of experimental data on particle deposition by Liu and 
Agarwal (O, Ref. 48) and by Wells and Chamberlin {-k, Ref. 47), with the results of 
our Montecarlo simulations (thick lines). The lower line corresponds to simulations 
without taking into account the effect of the Saffman lift. The thin line is the result 
of simulations from a one-fluid model (Ref. 39). 

grangian velocity defines a Lagrangian model obeying the Thomson-87 well-mixed 
condition, with uniform concentration PDF given by the Eulerian one-point PDF 
for the random field. 

Depending on the circumstances, a random field based approach to Lagrangian 
modelling may be advantageous. In the compressible case, knowledge of the one- 
point Eulerian PDF for the random field is sufficient to determine the coefficients 
of the associated Lagrangian model. The Thomson-87 approach, instead, requires 
knowledge of the particular Lagrangian PDF (indicated in j3] with g^) which orig- 
inates from an initial concentration profile equal to the instantaneous local fluid 
density, and which does not necessarily coincide with pe- Solid particle transport is 
an example in which implementation of the Thomson-87 approach is not straight- 
forward, unless ad-hoc hypotheses are made on the Lagrangian statistics. 

A second advantage of this approach concerns the non-uniqueness problem: 
knowledge of the two-point Eulerian correlations completely fixes the form of the La- 
grangian model, which is of interest for turbulent flows in complex geometry, where 
it is not clear which model satisfying the well-mixed condition, should be choosen. 
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The relation of some of the non-unique terms with hehcity [211 and rotation [Hj is 
confirmed, and additional terms associated with strain have been identified. Similar 
approaches, in which DNS informations on the two-point Eulerian correlations are 
used to determine the form of Lagrangian models, have been recently adopted in [2Z] . 
Alternative formulations for the treatment of the non- uniqueness exist, in which the 
Lagrangian acceleration is modelled by a Langeving equation, on the same footing 
of the Lagrangian velocity (second order models P). In these models, however, the 
non-uniqueness problem is only displaced to the higher order acceleration. 

A third advantage of the random field approach concerns situations in which it 
is difficult to characterize an inertial range, and in which concepts like the constant 
Co cease to be meaningful (e.g. in the buffer region of a turbulent boundary layer). 
Comparing our approach with the Thomson-87 technique, the main difference is, to 
lowest order in the SO (3) expansion, the form of the noise and the parameter a taking 
the place of Cq. Both noise expressions require knowledge of quantities estimated 
from large scale features of the flow: the viscous dissipation e and the Eulerian 
correlation time te, whose relative dependence (as the one between a and Co) is 
not an intrinsic characteristic of the models. In our approach, however, a precise 
relation can be obtained between the parameter a and the ratio of the Lagrangian to 
the Eulerian correlation time t^/te, which is valid also when the Reynolds number 
is low. Using for this ratio the value obtained in [10], we obtain a ^ 0.65. 

We have tested our model, with isotropic noise and Gaussian statistics, to study 
solid particle transport both in homogenous isotropic turbulence and in channel flow 
geometry. 

In homogeneous isotropic turbulence, the correct renormalization for the cor- 
relation time for the fluid velocity along the solid particle trajectories have been 
obtained without resorting to ad-hoc parametrizations. The form of Eqn. ()101|) de- 
scends directly from the random field and Markovianization along trajectories [see 
Eqns. (jHl) and ()65|)]. This illustrates the importance of the parameter a in providing 
the most simple characterization of space correlations in the turbulent fiow. It is 
important to stress that, had we not taken its contribution into account, Eqn. (jlUl|) 
would have been unable to reproduce the anisotropy of the time correlations. 



44 



An interesting aspect we have observed is satisfaction of the ergodic property in 
sohd particle transport by homogeneous turbulence, this, despite compressibility of 
the solid particle flow. It is not clear whether this is an artifact of the model; in any 
case, it is a non-trivial effect since the ergodic property can be shown to be violated 
by very simple compressible flows [see Eqns. ()7HI7H|1 and discussion therein]. 

In channel flow geometry, we have found good agreement with experimental data 
on particle deposition |Tfl I48j . and partial agreement with numerical data on near 
wall accumulation j^llHl- We stress that these results have been obtained without 
any parameter fitting, apart of the choice a = 0.65 inferred from |l0] . 

Clearly, a Gaussian model with isotropic noise cannot account for the effect of 
coherent structures and intermittency, which are an important feature of turbu- 
lence in channel flows. Imposing the appropriate form for the one-point PDF and 
going to higher orders in the S0(3) expansion allows consideration of these effect. 
Preliminary analysis suggests that inclusion in the model of non-Gaussianity, noise 
anisotropy and non-unique terms, strongly affects particle deposition and transport 
in wall turbulence. Modelling the structure of the turbulent correlations, based on 
empirical considerations, appears to lead to models that perform worse, compared 
to the data, than the simple isotropic Gaussian model, a situation similar to that 
observed in 49 . This suggests that careful consideration of the structure of the 
turbulent correlation, based on DNS data, may be necessary; this will be part of a 
different publication. 

Some issues remain to be clarified as regards the definition of a random field 
purely in terms of its local properties. The global extension provided by Eqn. (fT3jl . 
in which the random field is assumed "Markovian" along rectilinear cuts in space- 
time is only one of the possibilities. This choice produces effects on the form of 
the trajectories, which can be accounted for only in the non-Markovian approach 
described in Section VIII. (Markovianization along trajectories corresponds to con- 
sidering only local properties of the random field). An open question remains which 
global structure of a random field would lead, for fixed local structure, to trans- 
port properties which are approximated best by a Markovian Lagrangian model. 
This, beside understanding whether the assumption in Eqn. ()15|). which leads in 
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the Gaussian case to exponential scaling of the random field correlations, fits the 
turbulent structure in appropriate way. For the choice provided by Eqn. (fT^ . and 
for homogeneous isotropic conditions and Gaussian statistics, the error in the ratio 
Tl/te corresponding to a ^ 0.65 appears to be of the order of 15% in defect. 

Related to this issue is the fact that consideration of long-lived coherent struc- 
tures, corresponding to small values of tl/te and to large errors in the Markovian 
approximation, is probably out of the range of applicability of our model. An al- 
ternative strategy, which would allow taking into account long-lived coherent struc- 
tures, is the non-Markovian approach described in Section VIII. The noise and drift 
terms, however, would have to be rederived including the condition at the emission 
point following Eqns. dHZj) and (HHOl). 

Extension of the present approach beyond one-point statistics is possible in prin- 
ciple, but is limited by the unphysical scaling of the random field structure function 
at small separations. Only concentration fiuctuations on the scale of a correlation 
length could then be taken into account. 
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Appendix A: non-Gaussian case 

A symmetric one- dimensional distribution with unitary variance and kurtosis 
larger than three can be modelled by means of a bi-Gaussian: 



Pix) 



a 



— exp(-x ) + 



T12 



r exp - 

(27r/3)5 2(3' 



where 



a 



4{x^) - 12 



and P = -{x^) — 1; 
3 



A{x^) - 9 

parameterize the strength of the kurtosis (x^). From here, we can obtain the ex- 
pression for an isotropic non-Gaussian velocity distribution: 



PE 



aexpf ^ 



1 



a 

-exp( 



u 



2(3u. 



(A3) 
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and for an anisotropic distribution, in which one of the velocity components, in the 
diagonal reference frame for the Reynolds tensor, is non-Gaussian: 



where: 



Pe = api + (1 - a)p2 = Px[apyi + (1 - a)py2]pz (AA) 



Px = -, — i—T exp(--^), pz = - — T— T exp( 



1 

with i?P = Rf^ = (3R^'^, and the hat indicating the diagonal reference frame. 

Let us calculate explicitly the drift terms in the case of Eqns. (A4-A5). Substi- 
tuting into Eqn. (fT^ . we find immediately that A is given by the superposition of 
the contributions from each of the Gaussians pi and p2'- 

A;dx^ = -^^{dwW)[ap,S^^, + (1 - a)p,S],y (AQ) 

(notice the absence of the hats; it is not necessary here to work in the diagonal 
reference frame) . The contribution from the $ and terms has a more complicated 
form. Let us take the laboratory frame with the inhomogeneity direction along 
(the usual channel flow geometry in which is the mean flow direction). We use 
the ansatz: 

$; = a^;, + (1 - a)^;, + a$;, aI>; = P.s^ (a?) 

where $i and $2 give the form of $ in the case pE = Pi and pe = P2- Substituting 
into Eqn. ()2Up and using Eqns. (A4-A5), we obtain: 

du^F^ = -5ld2a{pi - P2) (AS) 

leading to the result in the laboratory reference frame: 

A$; = -^l^\^d2a / d^ e-"^/^ (A9) 

(27r)2 Ju2/y/W 

where f2*2 is the rotation matrix defined through = VL^ m\ i.e. fi* ■ = e* ■ e^-. 
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Analogous procedure is followed to obtain the \1/ term. From Eqn. ()21|). in 
analogy with Eqn. ((2S1), write: 



= a# + (1 - a)ij' + Alp' Alp' = S^G (AlO) 

where 2duiip^ = —A(^\ = —F2, and decompose \E'j* in analogous way. From Eqn. 
(A9), we obtain immediately: 



G = ^^n\d2a du dx e-^ /2 ^^^^^ 

2(27r)2 i_oo JulJ^ 



and, substituting again into Eqn. [in real space: ^ = S^dukip'' — 9„j '?/'*], we 

find: 

= S}d{,2G - Q\Q/di,iG {A12) 

Assuming that the statistics of dw be independent of the velocity, as we have done 
in section II, has the consequence that all of the non-Gaussianity is contained in 
the drift. The small scale structure of the correlations, associated with dw remain 
therefore Gaussian. This breaks, for large values of the kurtosis, the direct relation- 
ship between the drift coefficients and the correlation lengths It is easy to see 
what happens looking at Eqns. (A4) and (A6). Whenever the value of u goes above 
ut, the slowly decaying becomes dominant in Eqn. (A6) and leads to a reduced 
decay rate for the fiuctuation that can thus slowly grow to produce a burst. We 
thus have a hierarchy of time-scales (focus for simplicity on variations along time): 

Te — > Time scale for background {u ut) 

Pte — > Time scale for bursts {u ~ I3^ut) 

P'^te — > Spacing between bursts 

This difference between the time-scale for bursts and background fiuctuations is not 
physically meaningful in general. This has the effect of overshooting the burst contri- 
bution to turbulent dispersion, which is estimated by the product of the probability 
j3~^ of a burst, its time scale and the square of its velocity scale: 

X pTE X pul = pU^TE. (A13) 
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(In comparison, the background contribution is u'^te, and, if the time-scales for burst 
and background were the same, the background and burst contributions would be 
identical) . 

Appendix B: Spherical tensors and the SO(3) technique 

A symmetric two-index tensor function can be decomposed in spherical tensors 
in the form 

5*^Fj(x), d'd^Yj{^), x'x^Yj{^), {x'd^ +x^d')Yj{^), 

x„(e^""^a;* + e*""^a;^)5^F7(x) and x„(e-'"'"a* -|- e*""^a^ )9„Fj(x) (51) 
where Yj{x.) is a J-order polynomial in the components Xf. 

Yj{^)=y'^''-'-'x,,x,,...Xi, (52) 

and ?/*i*2 - v jg traceless with respect to any pair of indices ^H]- In consequence of 
this, the spherical tensors in Eqn. (Bl) will be polynomials of order L = J, J — 2, 
J + 2, J, J +1 and J — 1 respectively. In the case of the noise tensor (Aw'^Aw^), we 
have the additional symmetry with respect to spatial inversion, which imposes the 
condition that L be even. This implies J even for the first four spherical tensors and 
J odd for the last two. Limiting the analysis to J < 2, we notice immediately that 
the last spherical tensor in Eqn. (Bl) disappears. Similarly the J = 1 contribution 
from x„(e-''"™'x* + e*"''"x-^)9mYj(x) is absent due to incompressibility; writing Yi{x.) = 

9,x„(e^"™x^ + e*"'"x^>„Fi(x) = 5e^"™x„i/„ = 

which imposes ym = 0. We are thus left only with the J = and J = 2 contributions. 
From Eqn. ()39p. the J = and J = 2 contributions to S'-' will have the form: 

urBi^i^) = a\^\6'^ + a^] (53) 

UtBl^ix) = - — -XiXmS^^ + 2c'™|x|(9j9,X;Xm + -. — -;rXiXmX^X^ 

|x| jxl'^ 

+ —^{x^d' + x'&)xiXm (54) 
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Applying the incompressibility condition diB"^^ = leads to the equations: 

« = -! 

3c/'™ - 46'™ - e'™ = 

Substituting the solution to Eqn. (B5) into Eqns. (B3-B4) leads to Eqn. (jlU)). 

We give next the expressions for the spherical tensors contributing to an anti- 
symmetric two-index tensor: 

e'^'^XkYj e'^^dkYj and {x'&-x^d')Yj {B6) 

In the case of the antisymmetric part of the correlation Ca, we have the additional 
property of antisymmetry with respect to spatial inversion, which implies that J be 
even for the first two and odd for the last. 

In the case of a vector field, an analogous decomposition can be obtained in 
terms of spherical vectors in the form: 

x'Yj{^), d'Yj{^) and e'^%dkYj{^) (57) 

If the vector field does not have an axial component, only the first two spherical 
vectors can contribute. If we have axial symmetry, identified by a direction u, 
the tensors y*i'2---*' entering Eqn. (B2) will be zero trace symmetrized products of 
components and of the identity matrix 6^^ . The first spherical vectors x'^Yj{x.) 
and 9*Yj(x) are respectively: 

X, (u-x)x, ((u-x)2 - -luplx^x (S8) 

3 

and 

0, u, (u-x)u-i|upx {B9) 

3 

Appendix C: Conditional random field statistics 

We want to calculate conditional velocity moments in the form 

(m*(xo, to)M^(xo, to)---|u(Xi, ti), U(X2, ^2), •••) (CI) 
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Let us indicate, for 1 — 0,1, ...n, U; = u(x/,i/), and = (U^Um), assuming for 
simplicity a symmetric correlation tensor. For a Gaussian random field, the velocity 
correlations at points {t;,Xi} are obtained from the generating function 

p({77,}) = exp ( - - ^ r/, • C,^ • 77^) (C2) 

l,m=Q 

Let us introduce the marginal PDF 

1 ^ 

p({U,,i = l,...n})=AAexp(-- J] U,.D,^-U„) (C3) 



l,m=l 



where jV is the normalization and D;^ is the inverse of the restriction of Cim to 
l,m — 1, ...n. We shall indicate in the following this restriction with a prime: 

n 

^' ^ J] , {UJ} = {U,, / = 1, ...n}, etc. {CA) 

Im 1,171=1 

The generating function for Uq conditioned to U^, / = 1, ...n is obtained by inverse 
Fourier transforming p{{r}i}) with respect to {77^} in {U^} and dividing by the 
marginal PDF p({U;}): 

PMm) = -^^^ j exp ( - lY^m . U, 

-^^^l ■ ^Irn ■ »7m " Y^Vo ' " ^ - ^77o • Cqo ' T/q) (C5) 

im / 

From here we can calculate the conditional moments in Eqn. (CI). We calculate 
first the mean velocity in {to,xo} given velocities U/ in {t/,x;} / 1, ...n: 

1 ap(r/o|{U;}) 



(UolW}) = - 



p(0|{Ua) air7o 



rjo=0 



p({u;})p(o|{u;}) 



/ i 

Carrying out the Gaussian integrals, we obtain the result: 

(Uo|{U;}) = Y^C^i ■ • U„ (C7) 

im 

The calculation of the second conditional moment is analogous and the result is: 

1 d^p{'n,\{\j\}) 



(UoUoKu;}) 



T7o=0 



= Coo - Y^Coi ■ T>i^ ■ C^o + (Uo|{U;})(Uo|{U;}) (C8) 



im 
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